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FOREWORD 


The work reported herein was performed by the Grumman Aerospace Corpor- 
ation under the NASA/Langley Master Agreement and Contract No. NAS 1-10635 
for the Development and Implementation of Space Shuttle Structural Dynamics 
Modeling Technology. The Work Statement of Task Order No. 17, "Development 
of an Analytical Program to Analyze Reusable Surface Insulation for Shuttle", 
authorized and specified the tasks to be performed in this study. The period 
of performance was for 15 months starting in June of 1973. 

The overall supervision of programs under the Master Agreement is pro- 
vided by Mr. E.F. Baird, Master Agreement Program Manager. The Task Order 
No. 17 Project Manager was Dr. I. U. Ojalvo. Many individuals at Grumman 
contributed to the work reported here. However, the authors wish to speci- 
fically acknowledge the efforts of Mrs. Patricia Ogilvie for greatly assist- 
ing in the development of the associated computer program. 
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ABSTRACT 


This report describes an efficient iterative procedure for the vibration and 
modal stress analysis of Reusable Surface Insulation (RSI) of multi-tiled Space 
Shuttle panels. The method, which is quite general, is rapidly convergent and 
highly useful for this application.' A user-oriented computer program based upon 
this procedure and titled RESIST (REusable Surface Insulation Stresses) has been 
prepared for the analysis of compact, widely spaced, stringer- stiffened panels. 
RESIST, which uses finite element methods, obtains three diemensional tile 
stresses in the isolator, arrestor (if any), and RSI materials. Two-dimensional 
stresses are obtained in the tile coating and the stringer- stiffened primary- 
structure plate. A special feature of the program is that all the usual detailed 
finite element grid data is generated internally from a minimum of input data. 

The program can accommodate tile idealizations with up to 850 nodes (2550 degrees- 
of-freedom) and primary structure idealizations with a maximum of 10,000 degrees- 
of-freedom. The primary structure vibration capability is achieved through the 
development of a new rapid eigenvalue program named ALARM (Automatic LArge Reduc- 
tion of Matrices to tridiagonal form) . 
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I. INTRODUCTION 


The integrity of a reusable space shuttle system is strongly dependent 
upon protecting the orbiting vehicle from reentry heating in a manner which 
does not require significant refurbishment between missions. The thermal 
protection system (TPS) selected to fulfill this need is a non -structural re- 
usable surface insulation (RSI) material shaped into individual tiles which 
cover almost all of the orbiter surface area. Most of the tiles are square in 
planform and measure 15*3 or 20.4 cm (6 or 8 inches) on a side, with thicknesses 
which vary roughly between 1.3 and J.6 cm (0.5 and 3 inches). Since loss of a 
single tile could be catastrophic, RSI tile stresses must be accurately deter- 
mined for the anticipated static, dynamic and thermal environments. 

Recognizing that a system' s natural vibration modes may serve as basic 
building blocks for predicting the response to acoustic excitation, this re- 
port describes a new iterative procedure for accurately determining tile 
stresses associated with typical shuttle-panel lower frequency modes. The re- 
sults of this work must be combined with a dynamic response program, to obtain 
realistic forced response results. 

A user-oriented computer program based upon the present method of analysis 
was developed. The program, which is titled RESIST (for REusable Surface 
Insulation STresses), is capable of computing undamped vibration mode shapes 
and frequencies of elastically supported multi-tiled panels and determining 
associated normalized RSI stresses. Typical numerical results from this com- 
puter program are presented. A user’s manual to facilitate its implementation 
is presented as Appendix B. 

Because of the complex geometry, nonuniform anisotropic material properties, 
and detailed three-dimensional stress states, the TPS was idealized by finite 
element assemblages with up to 2500 degrees of freedom per tile. Since a 
number of tiles affixed to a given structural panel will, in general, interact 
with one another, application of the standard direct stiffness method would 
require simultaneous equation systems involving excessive numbers of unknowns. 

The present iteration scheme overcomes this problem by treating each tile 
separately. An important byproduct of this approach is that it avoids condi- 
tioning problems associated with combining low-stiffness tile isolation ele- 
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merits (Es 50 psi) with high-stiffness primary structure elements 

S’ 

(E~10 x 10° psi). Typical results from the associated computer program reveal 
a rapid rate of convergence. In addition, a related effort for obtaining tile- 
stresses associated with statically loaded and heated shuttle panels is pre- 
sented in Reference 1. 
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II. TECHNICAL APPROACH 


A. STRUCTURAL CONFIGURATION 

The design configuration for which an analysis procedure is presented is 
shown in Figure 1. It consists of a stringer-stiffened flat rectangular panel 
which supports a nonstructural thermal protection system (TPS). 

The TPS is composed of a series of rigidized (RSl) tiles. The tiles are 
undercut on all four sides to accommodate "filler-strips". The purpose of the 
nonrigidized filler-strip insulation is to prevent severe heat penetration 
through the gaps between adjacent tiles. 

The RSI material is not bonded directly to the primary structure, but to 
a thin, stiff, strain arrestor plate first, and then, in turn, to a soft strain 
isolator material. The function of both these items is to help isolate the 
primary-structure strains from the low-strength RSI tiles, for the wide range 
of loading environments which the vehicle must sustain. 

B. GENERAL SOLUTION PROCEDURE 

Because of the detailed complexity of the structural configuration and the 
dependence of material properties upon temperature, an analysis technique 
based upon finite element methods was selected. However, direct application 
of the standard direct stiffness finite element procedure, to obtain accurate 
three-dimensional tile stresses, would require equation systems involving’ ex- 
cessive numbers of unknowns. The reason for this is that many tiles affixed 
to a given structural panel may interact. Since each one is a three-dimensional 
body requiring a detailed structural idealization, their simultaneous conside- 
ration requires the solution of large systems of equations. To overcome this 
problem a rapidly convergent iteration scheme, which treats each tile separate- 
ly, was developed. The logic flow for an efficient computer program which 
employs this procedure is presented in Figure 2. The basis for the method is 
that the TPS is nonstructural but its stress levels, which are critical, mutt 
be computed. Thus, it becomes possible to neglect the ’stiffness of the TPS 
initially, but not its mass, to determine approximate primary structure modal 
deflections associated with the vibrations of the overall system. In these 
initial calculations, the tiles are assumed rigid in shear and thickness- 
stretching for the purpose of computing their kinetic energies. 
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Figure 2. Logic Diagram for Solution Procedure 
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An interative procedure is then performed where, for each step, the primary 
structure oscillatory modal deflections are imposed individually upon each, tile 
at the tile/primary-structure interface, and the tile deflections, and interface 
boundary loads are obtained. The frequency is then updated by computing a 
Rayleigh Quotient, using the latest non-rigid tile displacements in addition to 
the corresponding primary structure displacements. The individual tile 
boundary loads obtained are then assembled and their reactions applied to the 
primary structure. New primary structure deflections are obtained ‘and compared 
to the previous set. This process is repeated until convergence is obtained. 

By this method each tile is temporarily assumed uncoupled from all others. Al- 
though this is not strictly true', the coupling involved is sufficiently weak 
so as to ensure accurate approximate results and rapid convergence as well> . 

It should be noted that an important byproduct of the present tile-by-tile 
solution procedure is that, unlike the direct stiffness method, it avoids 
possible numerical precision problems associated with directly combining a low 
finite element stiffness (E = 50 psi) with a high primary structure stiffness 
element (E = 10? psi). ■ 

C. STRUCTURAL IDEALIZATION 

A pictorial representation of the four finite element types incorporated 
in the RESIST computer program is presented in Appendix C. A brief description 
of these elements and how each type is used in the overall structural ideali- 
zation is presented below, y • - . . 

Primary Structure 

The primary structure stiffness, idealization is based upon two finite ele-? 
ments contained in the Grumman program library described in Ref 2. The panel 
surface is represented by flat rectangular membrane-bending elements. The ele- 
ments possess three deflections and two rotations at each node, for a total of 
20 independent nodal degrees of freedom per element. The stiffeners are re- 
presented by beam elements which bend, stretch and twist. These elements pos- 
sess 6 degrees of freedom per node, for a total of 12 degrees of freedom per 
element. •: ■ ' ; . 

In general, the stiffeners do not attach to the plate nodes at the beam 
centroids. The stringer attachment points actually used coincide with the plate 
node which is closest to the centroid of each beam element. The assumption 
associated with this connection is that the attachment point and beam centroid 



move as two points on the same rigid body. Furthermore, the principal axes of 
the beam cross section may make an arbitrary angle with the plate surface. Re- 
fer to Appendix C for further details on the manner of attachment of beam ele- 
ments to plate nodes. 

The translational mass properties of each rectangular plate element are 
divided equally and lumped at their middle surface node points. The stiffener 
translation mass properties of each beam element are also divided equally 
and Tumped at their two centroidal nodes. The rotational inertias of the 
primary structure beam and plate elements are ignored. The assembled struc- 
ture is rectangular in planform and is assumed to be statically supported 
with uniform, but different, arbitrary boundary conditions along each edge. 

The above idealizations and assumptions appear consistent with the degree 
of detail required for obtaining accurate dynamic responses of typical aero- 
space stiffened panel construction. 

Thermal Protection System (TPS) Tiles 

The basic stiffness element used for the RSI tiles and strain isolator 
is the anisotropic hexahedron, as defined in References 3 and 4. The version 
of this element which is being used contains 3 deflection degrees of freedom 
per node and 8 nodes per element (the more general version which is available 
permits up to 20 nodes per element). A typical undercut tile with element 
and node numbering is shown in Figures 3 & 4. The top layer of each tile also 
has a coating material which is idealized by rectangular membrane elements. 

The bottom two layers of each tile comprise the strain isolator and strain ar- 
restor, respectively. The undercut regions of each tile are filled with non- 
rigidized insulation material. This contributes a small amount of mass to the 
structure but negligible stiffness and negligible effect upon the stresses 
within the rigidized tile material. Thus, for purposes of computing the 
tile deflections and stresses, the undercut regions are treated as empty voids. 

Mass Matrix Calculations 

Two different mass matrices for the TPS must be computed. The first mass 
matrix required for the TPS "HEX" element and used in Steps 4 and 5 of the logic 
flow diagram (Figure 2),' is based upon lumping l/8 of each element's mass at 
each of its 8 corner nodes for x, y and z translational motion. The second type 
mass matrix, required in Step 2 of Figure 2, is necessary for computing the ap- 
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Figure 3. Typical Finite Element Idealization of Shuttle RSI Tile 
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Figure 4. Top View of Lowest Surface of Deck 
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proximate mode shapes and frequencies (Step 3 of Figure 2 ) based upon certain 
assumptions. These assumptions are that the TPS behaves as a series of rigid- 
ly attached masses which contribute inertia, but afford no structural stiffness 
to the primary structure. Analytical details associated with the latter pro- 
cedure are described as follows. 

Let u, v and w be the displacements of a generic TPS point above a 

corresponding point on the primary structure plate (p) with displacements 

u , v and w . Then under the assumption of a rigid-in_shear and thickness- 
P P P 

stretch TPS we obtain 



z 


P 

Sx 


v = v 

P 


3w 


z 


ay 


w = w 

p - 

(i.e., normals to the plate middle surface remain normal). 

For a section of TPS above a typical plate element, the kinetic energy 
is given by 
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where A is the planform area surrounding a node. 
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Defining the integrals in the kinetic energy equation as follows (refer 
to Figure 5): 



and computing 

d_ 

at 


then yields for the TPS mass matrix, per plate node 




2-9 



ARRESTOR (p 2 ) 


Z = Z 3 - z 2 + t a 



Z = Z 2 = Z V + T F 


Z = Z 1 



Z = 0 


Figure 5. TPS Z-Coordinates for Computation of Rigid Mass Matrix 



Automatic Grid Generators 


To reduce the quantity of input data significantly, two automatic mesh 
generating schemes were developed. The first procedure generates the primary 
structure geometry, node and element numbering system, material properties, 
and boundary conditions; while the second does the exact same thing for the 
TPS. Each requires a minimum of information, such as overall dimensions 
stringer pitch, material properties (as a function of temperature) and TPS 
material temperatures and element properties from which data may be inter- 
polated (see Appendix B for details). 

D. APPROXIMATE FREQUENCIES AND MODE SHAPES 

This phase employs the primary- structure mass and stiffness matrices (the 

T 

latter being in lower triangular decomposed form L, where K = LL ) as well as 
the rigid-tile approximate mass matrix, just described, to obtain approximate 
mode shapes and frequencies. Since the primary structure idealizations are 
quite refined, the number of degrees of freedom will be in the thousands 
with 10,000 as an upper limit. Thus, an efficient and reliable eigenvalue 
reduction scheme is required. 

The reduction procedure employed was suggested by Crandall^ and is 

based upon an algorithm for accommodating the full matrix due to Lanczos^. 

However, the basic scheme is numerically unstable. - The essential improvements 

(7) 

necessary to correct this weakness were made by Ojalvo and Newman , who 
were the first to develop a successful reduction scheme for large scale pro- 

/O \ 

blems. The work was further streamlined by Newman and Pipano v .through the 
use of banded matrix packing techniques and the computation of convenient 
error bounds on frequency squared. 

The present program, which is a synthesis of ...these works, is called ALARM 

Y 

(Automatic LArge Reduction of Matrices). It is designed to operate efficiently 
through use of the Grumman Comprehensive Matrix Algebra Procedures (COMAP) in- 
terface programs developed by C. Wilkie and F. Nolan, 
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A mathematical statement of the original vibration problem considered 
is 

[K] [ 0 ] = * [M] [ 0 } ( 1 ) 

where n x n stiffness matrices [M] and [K], which are generally highly banded 
and sparse, may be singular but non-negative, with the only restirction be- 
ing n ^ 10,000. The associated reduced eigenvalue problem is symmetric, tri- 

( 7 ) 

diagonal and of the. form 



2 

where m is typically several orders of magnitude smaller than n, and \ = V tu . 

The algorithm used for obtaining the a. and B. of Eq. (2), as will be 

(9) 1 1 

shown, is based upon the power' method . • Unlike the power method, how- 
ever, in which each vector is discarded as a new one is generated, the pre- 
sent procedure retains each vector after orthogonalizing and normalizing 
it to all previous vectors such that ]_ V. J fV.} = where 6 . . is the 

Kronecker delta. Thus, the involvement of the original system's lower modes 
is strengthened. A Rayleigh- Ritz approximation, using the' n x m matrix 
[V] whose columns are the V^, follows from the substitution of 

[L] T [ 0 ] = [V] [y] ( 3 ) 

into Eq. (l) where 

[K] = [L] [L] T (k) 


2-12 



Premultiplication of the resulting equation by K [V] [L] then yields 
Eq. (2). This result is now solved (in core) by Sturm- Sequencing for the 
X's and power iterations (as described on page 622 of Reference 9) for the 
eigenvectors. 

(7) m 

Experience has shown v that over half the number of modes, ^ ? obtained 

by this method are "exact" with a sharp drop-off in accuracy somewhere around 

p 

thus, the program usually sets m to 2q+l, where q is the number of desired 
3 /g\ 

mcdes. Newman- Pi pano error bounds v ' are useful indicators of where a loss 
of accuracy occurs. The computational steps for obtaining the approximate 
frequencies and mode shapes are presented in the Appendix A. 

E. FLEXIBLE TILE SOLUTIONS 

The tiles are originally treated as nonstructural mass items. However, 
since their modal stress states are desired, it is necessary to eventaully 
account for their flexibility. This is achieved by imposing the primary- 
structure/TPS interface deflections {figj^sin out, to each tile and computing 
the associated tile deflections, ft,.}, sin out, to obtain tile stresses and 

t\ 1 

reactions, .{-Pgj sin out-.. As described earlier, this, is repeatedly done in ~ 
an iterative manner until, convergence is obtained. 

Referring to Figure 6, the associated steady-state partitioned matrix 
equations governing each tile are 



where use has been made of the fact that the tile mass matrix is diagonal. 
This equation is then' decomposed into 

lK m M M ] [6 A }. = .[K^] (6 B } i (6) 

and 

''BA-* ^'A^i " ^ K BB _U) ^B^ ^ 6 B^i = ^" P B^i (7) 
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Prior to solution of Eq. (6) for {6^}^, it is noted that [K^ -ui 
is not positive definite if any of the eigenvalues, \^, of 

X A ^ £‘a> = E M AA^I {6 a 1 


satisfy the equation 


X A* 


Cl) 


(8) 


where o> is the current (approximate) frequency value. 

Since the efficient solution routine^) that is used for large banded 
symmetric equations, such as Eq. (2), requires positive definiteness, Eq. (6) 


is premultiplied by (K 


AA 


-CJU 


M^) to insure this property. Thus, the equa- 


tion treated to obtain £6 is 


[K aa A.j Wh * [K m {-K m 6^ 


(9) 


AA “AA" 

Following this, Eq. ( 7 ) is used to compute [-P^]^* Once convergence of 

the primary structure deflections has been achieved, the corresponding tile modal 
stresses are computed from the hexahedron stress recovery equations ( 3 ). 

F. FREQUENCY UPDATE 


Since the modes and frequencies obtained originally ignored the tile 
stiffness and true displacements, the resulting frequencies are only approx- 
imate. To correct these, the frequencies are updated through use of a Ray- 
leigh Quotient calculation in which the latest tile and primary structure 
deflections are used in the equation 


2 5 l6tj i m ( H + L6ps J [v] M 

U) = 7 — • (10 

i LTJi W M,* L6 pe J [M p< ] { 6 ps }_ 

where the subscripts "T" and "ps" denote tile and primary structure, respec- 
tively. Thus, (6 } are the latest primary structure deflections, which 

ps 
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originally are the {0} of Eq. (l), [K J and [M ] are the banded stiffness 

ps ps 

and mass matrices, respectively, of the primary structure and tile. Thus, 
the [Krj,] and [M^J matrices of Eq. (10) are related to the stiffness and mass 
matrices of Eq. (5) as follows: 


k aa' 

k ab 


^B 


[k t ] = I 


M 


AA 0 
0 


G. PRIMARY- STRUCTURE DEFLECTION UPDATE 

The original primary- structure modal deflections are approximate, should 
be checked and, if necessary, corrected. Employing the latest compatible set 
of tile and primary- structure deflections, the frequencies are recomputed, to 
give by the method described in the previous subsection. The 

steady oscillating tile boundary reactions sin u> t are collected for all 

the tiles and applied to the primary- structure, to yield new modal deflec- 
tions, as follows: 


[K -«J 2 Ml f 6 } = {-PJ (11) 

L ps R ps - ps J 1 B i 

where f-P-r,} sin uo D t is the assembled load reaction of all the tiles acting 

.D -K 

upon the primary- structure. 

( 2 ) 

Once again, since the solution routine that is used for large, 

banded symmetric equations requires the positive-definiteness property, 

2 

Eq. (ll) is premultiplied by -(ju p M ] prior to solution for the 

{6 ps 3, i.e., the actual equation solved for {6p g } is: 


O C P 

[K M ] {6 } = [K -u)„ M ] {-P ] 

L ps R ps c ps J ps R ps 1 ps 


( 12 ) 


For consistency, the deflections { 6 p S ) are then normalized such that 


L6 J [6 } = 1 

ps J l ps J 
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prior to comparison to the previous set, to establish convergence. The con- 
vergence test employed consists in satisfying 


max element of 


6 -6 
Ps Ps 

present previous 


£ e 


largest element of 6 


ps 

present 


where g is an empirically determined positive input quantity. 
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Ill . NUMERICAL RESULTS 

A. TILED CONFIGURATIONS 

To demonstrate the numerical convergence properties of the proposed 
iteration scheme, two examples are presented here. The material properties 
and loading conditions used, while not precise in a specific design sense, 

are representative of the problem parameters the method is intended to 
accommodate. 

Figure 7 shows the configuration, finite element idealization, mater- 
ial properties, and boundary conditions considered for Example 1. The 
fundamental frequency, largest primary- structure normalized deflection 
and average stress components for certain critical RSI elements are pre- 
sented as a function of iteration number in Table 1. Besides revealing 
a rapid convergence rate, the results presented indicate a high level of 
accuracy for the first iterate in that the maximum frequency, deflection and o x 

errors are less than 2 %. The errors in a were higher for the first itera- 

z ■ * 

tion (approximately 20%), but these settled down to under h% in the second 
iteration,, 

Similar observations were made for the more realistic shuttle tile con- 
figuration of Example 2 shown in Figure 8. This problem consists of three 
6x6 inch tiles, each of which is two inches thick and has 1 x l/2 inch 
edge undercuts. The primary structure is a .Oil inch aluminum plate with 
offset stringers spaced half an inch apart. The two opposite shorter 
boundaries are clamped while the longer ends are free. 

The fundamental frequency and peak normalized deflection, as a function 
of iteration number, are shown in Table 2. The primary- structure mode- shape 
is plotted in Figure 9 for iterations 1, 2, and 4. The corresponding stresses 
for iterations 1 and b are shown in Figures 10 through 12. Running time for 
this realistic type problem on an IBM 370 averaged. 1 l/3 CPU minutes per tile 
plus a fixed time of 13 minutes per iteration. Thus, four iterations of the 
three tiled configuration required approximately 68 CPU minutes. 
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Figure 7. Example 1 : Sample Problem Used to Demonstrate Numerical Convergence of Iteration Scheme 



TABLE 1 


FUNDAMENTAL FREQUENCY, MAXIMUM DEFLECTION, AND CRITICAL 
TILE STRESSES VS. ITERATION NUMBER (Refer to Figure 7) 
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Figure 8. Example 2: Typical Configuration That Computer Program is Capable of Analyzing 



TABLE 2 


FUNDAMENTAL FREQUENCY AND MAXIMUM PRIMARY STRUCTURE 
DEFLECTION AS A FUNCTION OF ITERATION NUMBER FOR 

EXAMPLE 2 


ITERATION 

NUMBER 

FUNDAMENTAL 

FREQUENCY 

(Hz) 

NORMALIZED 
MAX PRIM STRUCT 
DEFLECTION 
(inches) 

0 

l46 

0.111 

1 

153 

0.108 

2 

153 

0.111 

3 

154 

0.110 

4 

153 

0.110 
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Figure 9. Example 2: Plate Structure Normal Deflections of First Mode 
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Figure 10. Example 2: Normalized First Mode Direct Stresses for Tiles 1 and 2 in Planer Directions as a Function of Spanwise Coordinate 






Figure 11. Example 2: Normalized First Mode Direct Stresses for Tiles 1 and 2 
in Planer Directions as a Function of Spanwise Coordinate 
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PLANE OF SYMMETRY 




Figure 12. Example 2: Normalized First Mode Direct Stresses for Tiles 1 and 2 as a Function of Chordwise Coordinate 
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B 


UHTILED PRIMARY STRUCTURE 


Figure 13a is representative of a typical stiffened untiled "shuttle panel" 
■undergoing dynamic testing at the Langley Research Center. Figure 13b shows the 
finite element idealization that the RESIST computer program generated for this 
panel. The numerical results obtained (see Table 3) reveal a tendency for the 
frequencies to occur in groups of 7. The geometric significance of this clustering 
is that the modes in each group are related to the physical presence of 7 stringers, 
which subdivide the panel into 6 bays. ■ 

The idealized stringers are assumed to be attached at only a single rivet line. 
However, in the configuration of Figure 13a, each stringer is attached by two rivet 
lines to form a closed section torque-tube. Thus, the idealization of Figure 13 b 
is inadequate for predicting the higher cross-stringer modes (i.e. m >0) when the 
stiffened plate-pitch, a, is not small compared to the unstiffened plate-pitch, b. 
Since a is of the order of b for realistic shuttle panels, the idealization capa- 
bility within RESIST should be extended to accommodate closed-section stringers 
with narrow spacing in addition to the present capability which treats widely 
spaced, or open-section, stringers. 
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IV. CONCLUSIONS AND RECOMMENDATIONS 


An iterative procedure for the vibration stress analysis of RSI multi- 
tiled shuttle panels has been developed. The method, which is quite general, 
is rapidly convergent and highly useful for this application. 

A user-oriented computer program based upon this procedure and titled 
RESIST has been coded. RESIST, which uses finite element methods, obtains 
three dimensional tile stresses in the isolator, arrestor (if any) and RSI 
materials. Two dimensional stresses are obtained in the tile coating and the 
stringer- stiffened primary structure plate. A special feature of the program 
is that all the usual detailed finite element grid data is generated internally 
from a minimum of input data. 

The program may be used in an iterative mode to obtain detailed results. 
.However, for parametric design studies, reasonably accurate results may be obtained 
by using only one iteration at a significant savings in running time. This is 
achieved by having the program analyze certain specific tiles of a multi- 
tiled panel only once, while ignoring less critically stressed tiles. 

At present the program can accommodate tile idealizations with up to 
850 nodes ( 2,550 degrees-of- freedom) and primary structure idealizations with 
a maximum of 10,000 degrees-of- freedom. In addition, the tile pattern must 
begin and end at a structural frame and the isolator material must be iso- 
tropic. Should such restrictions, or any similar ones, require alteration as 
the shuttle TPS design changes, it would appear logical to update RESIST. 

To enhance the usefulness of the present work, it should be extended to ac- 
commodate closely spaced, closed-cell stiffeners (which are more typical of Space 
Shuttle panel designs) and to perform acoustic response analyses as well. 
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V. NOMENCIATURE 


A 

[A] 


[G] 

[K] 

CV 

^ K ps ] 

[ K AA^ ^ k ab^’ 


[L] 

[M ] 

L ps J 


CHp] 

[M] 


[M] 

^- m aa^ ^b^ 
t- p B>' 



'V 

tv 


Plate area associated with a given primary- structure node 

-1 -T 

Matrix used in large eigenvalue algorithm (= [L]~ [M] [L]~ ) 

Matrix used in large eigenvalue algorithm 

Stiffness matrix used in large eigenvalue algorithm 

Tile stiffness matrix 

Primary structure stiffness matrix 

Matrix partitions of [K^,] 

Lower triangular decomposition of stiffness matrix 
Diagonal primary structure mass matrix 
Mass matrix associated with rigid tile 
Tile mass matrix 

Approximate tile mass matrix associated with a typical plate 
node 

Mass matrix used in large eigenvalue algorithm 
Diagonal partition matrices of tile mass matrix 

Reaction loads caused by all tiles acting upon primary 
structure 

Load acting upon i^* 1 tile caused by primary structure 

Approximate kinetic energy associated with a typical node 
Unnormalized and unrefined Lanczos vector 

Refined but unnormalized Lanczos vector 
Normalized and refined Lanczos vector 
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z. 

1 


f V 

ts t ) 

m 


n 

n 

q 

t 

U, V, w 



Thickness coordinates. Refer to Figure 5. 

Vector used in eigenvalue algorithm 
Vector used in eigenvalue algorithm 
Order of reduced eigenvalue problem 
Order of unreduced eigenvalue problem 
Number of dynamic degrees of freedom 
Number of desired approximate system modes 
Time 

Approximate displacement components of arbitrary point 
within tile 

Displacement components of primary structure 


y, z 
{y}> 

m 


y i 


or. , 8 . 

i’ l 




[6 A ]i> (Sgli 


f6 T 1i 


e 

X 


Spacial coordinates 

Eigenvectors of reduced eigenvalue problem 
Last element of {y^} 

Diagonal and superdiagonal elements of tridiagonal reduced 

eigenvalue matrix 

Primary structure deflections 

Deflections of i^ tile at node points not in contact (A) and 
in contact (B) with plate 

Deflections of i^* 1 tile 

Primary structure deflection' convergence parameter 
Eigenvalue of reduced problem 
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p 


Eigenvalue problem associated with "A" degrees- of- freedom of 
tile 

Mass density 


P n i 
f 0 } 

- ^ps ' 


Integrated mass density parameter 
Approximate primary structure mode shape 


in Frequency of large eigenvalue problem 

2 

l0 R Rayleigh Quotient 

( U Approximate rigid-tile frequency 
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APPENDIX A 


LARGE EIGENVALUE SOLUTION ALGORITM 


The algorithm employed to initially determine the approximate mode 
shapes and frequencies of the system^?*®) . j_ s outlined below: « 


1. Initialization: 


2. Generate (V Q ] : 


3. Compute 0 i ^ 


It. Test i s m+1: 
5. Compute [V^} = 


Set i = 0 and determine ro, the order of the reduced 
problem as follows: 

For n ^ 7 set m = n 

For n > 7 set m = n if 2q + 1 > n 

m = 2q + 1 if 2q + 1 < n and q > 3 

m = 7 if q £ 3 

"Reverse", transpose, and "reverse" [L] to obtain 
RTR 

[L] and save result for repeated use in Step 6 
(see discussion in Step 6). 

Use a random number generator for each element. 

Go to 6. 


iV, j 


- 1 - tv 

if no--continue; if yes--go to 15 


‘i-1 


{V_^}: This step normalizes the vectors such that 


LV, J (V = 1 

6. Solve for Cg^i and save it if i s 1. This is the solution vector of 

. [l] t {g.} = (V.} 

The usual Grumman routine used for solutions of 
this type equation is titled QBAC (see Reference 
10). However, it is not too efficient when there 
is only one right-hand side, as is our case. 
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A more efficient procedure is to 
manipulate [L] , which is an upper triangular 
matrix, into lower triangular form, and then to 
use QFOR to solve the equation 

[L] RTR [g i } R = {V.} R , where {g i ) R and { V ± } R corres- 
pond to (g^j and (V^j, respectively, each with the 

order of its elements reversed (e.g., the first 
element of {g^} is the last element {g^} and vice 

T RTR 

versa). The manipulation of [L] into [L] re- 
quires two matrix row reversals (COMAP routine 
REVERS*) and one transpose (TRAN*) operation in the 
order REVERS-TRAN -REVERS. The most time consuming 
routine of these three is TRAN, which is even 
slower than QFOR*. However, since it is only nec- 
essary to compute [L] once (see Step l), each 
solution for [g^I, by the present procedure, re- 

quires far less running time than by the direct 
method (e.g., l/30th the CPU time for a 3,000 DOF 
system) . 

7. Compute {b^} from {b^} = [M] (g }. This is actually accomplished by 

■ forming ( |_g^J [M]) which is considerably more 

efficient using the COMAP associated subroutines. 

8. Solve for (A V^} from [L] (A V^}=(b^} using QFOR. 

Note that [A] = [L] R [M] [L] Thus (A V^} is 
one power iteration of and so accentuates its 
lower modal content. 

9. Test i: if i = 0, set V^ = A and return to step 3. 

if i V 0, continue. 


*See Reference 10 for descriptions of these subroutines. 



10 . 

11 . 
• 12 . 


13. 

14 . 

15. 

16. 


Compute = |_V^J (A V^} and step up i by 1. 


For i = 2, compute 


This step computes the Rayleigh Quotient associated 
with the vector (V.}. 

- (1) l 

V 2 = (A V } -« 


For i s 3j compute 


(- X 1 )] 

v. 


Iav.1 

- O'. n 

[v., 

1 

‘ - B. o 

1 1 


l i-M 

i-1 

t "- 1 

1 i-2 


Test s > s: 


V i-2 


and set s = 1. 

This step theoretically^ orthogonalizes [V^}, 

However, any slight numerical inprecision tends to 
bias the resulting vestors unfavorably through the 
present iterative generation scheme. The suggested 
correction to this numerical stability problem is 
contained in the following two steps. 

if no, continue; if yes, set (V^} = 

return to 3. 


(s)| 


V. 

l 


and 


(s) 


Replace j 


by \ V 


“I - tv h“’ !){’,} 


for j = 1 


( s) | 


Continue correcting j ^ by repeating the above 

step for j = 2, 3, ...., 1-1. Step up s by 1 and 
return to 13 . 

"fchs 

Solve for the lowest 3 A (rounded to an integer) of the frequencies 


of Eq. 2 using Sturm sequencing 


(9) 


Solve for the eigenvectors associated with the highest "q" eigenvalues 

X, each normalized such that 

ly J Cy] = i 
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17. Compute and print m. = ' and f. = ^ 

11 1 2 tt 


|S y . m 

18. Compute — — j where y^ m is the (last) component of {y^}, 


X i * 


This term is an error bound 


(9) 


on 


<u i 


-1 


“i 


exact 


19. Compute [G] (y^j for i = 1, 2, , q where the columns of [G] are the 

vectors [g^}. These are the vibration mode shapes 

{0^} associated with the lower q eigenvalues. 


These are then normalized such that 
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USER'S MANUAL FOR 
RE*S*I*ST 

(STATIC AND DYNAMIC REUSABLE SURFACE INSULATION STRESS PROGRAM) 
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A. INTRODUCTION 


This Appendix describes the use of a finite element based structural com- 
puter program for determining the static response and natural vibrations of 
TPS protected shuttle panels. The program is titled "RESIST" for static and 
dynamic REusable Surface Insulation Stresses. The logic flor for RESIST is 
presented in Figure B-l. 

The basis for the method is that the TPS is nonstructural but its stress 
levels, which are critical, must be computed. Thus, it becomes possible to 
neglect the stiffness of the TPS initially, but not its mass in the vibration, 
to determine approximate primary structure deflections. 

An iterative procedure is then performed wheije, for each step, the primary 
structure deflections are imposed individually upon each tile at the tile/ 
primary- structure interface, and the tile deflections and interface boundary- 
loads are obtained. For the vibration option, the frequency is updated by com- 
puting a Rayleigh Quotient, using the latest non-rigid tile displacements in 
addition to the corresponding primary structure displacements. The individual 
tile boundary loads obtained are then assembled and their reactions are applied 
to the primary structure. New Primary structure deflections are obtained and 
compared to the previous set. This process is repeated until convergence is 
obtained. 

B. PROGRAM LIMITATIONS 

The usual assumptions for programs based upon the linear elastic finite 
element method are applicable to RESIST. However, to facilitate the prepara- 
tion of program input, a number of simplifications regarding the configura- 
tion and loadings have been made. Thus, the generation of a voluminous quan- 
tity, of finite element input data has been greatly reduced by inclusion of 
a series of data preprocessing subroutines within RESIST. The restrictions 
upon which these subroutines are based follows: 

1. Boundary conditions and edge loadings are assumed uniform along 
the four rectangular plate edges defined by x = 0, L and y = 0, 

L . 

y 

2. The primary structure plate temperature and properties are all 
uniform. 
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II 


MESH 


FORM K 

GENERATION 


& SET i = 1 



FIRST 

STATICS 

ITERATION 


COMPUTE TILE 
THERMAL & 
EXTERNAL MECH 
LOADS 


^ ST ? VB N 
(STATICS) 
(VIBRATION) 





COMPUTE 



RAYLEIGH 

4 — - ST ? VB \ 

FORM MASS 

QUOTIENT 

\ 

MATRICES 


VI (HEX) 


V (ALARM) 


YES/ GET APPROX 
RESULTS 



VES (must be ST prob) 


I NO (may be ST or DV) 


SOLN 

CONVERGED 
. ? 


i > i MAX 
? 


PRINT 

WARNING 

MESSAGE 


COMPUTE 

STRESSES 


COMPUTE 

STRESSES 


* APPLY FOAM/PRIMARY-STRUCTURE BOUNDARY DISPLACEMENT TO 
EACH TILE SPECIFIED, ONE AT A TIME, AND COMPUTE ALL NON 
BOUNDARY TILE DISPLACEMENTS AND STEADY STATE BOUNDARY 
LOADS. 

** FOR VB PROBLEMS USE LOADS FROM VI AND FREOUENCY FROM 
VII. FOR ST PROBLEMS WITH i = 1, USE LOADS FROM MODULE X. 
HOWEVER, IF i A 1 GET STATIC TILE LOADS FROM VI AND BOUN- 
DARY LOADS FROM X. 


PRINT 

OUTPUT 

& 

STOP 


Figure B-1 Flow Chart for RSI Stress Analysis Program "Resist' 
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3. The stringers are equally spaced with temperatures and properties 
which are all uniform. 

4. All tiles are geometrically identical as are their temperature 
distributions and uniform pressure loadings. 

5. The boundary conditions must be selected such that the primary 
structure is statically stable. 

The remaining limitations are primarily concerned with the program's 
' capacity and should be adhered to by the user. These limitations are as 
follows : 

6. Maximum number of nodes in a tile = 850. 

7. Maximum number of finite elements running in any one direction 
in a tile = 20. 

8. Ma xim um number of nodes in primary structure = 2500. 

9. Maximum number of primary structure nodes along x or y direction = 

1 , 000 . 

10. Maximum number of degrees of freedom in primary structure 

10.000 for vibration option. 

15.000 for statics option. 

11. Maximum number of natural mode shapes = 50. 

12. Maximum number of stringers =15. 

13. To avoid a singular stiffness matrix, I , and sin p must not 
both be zero for a given stringer. 

A violation of restrictions 6-13, inclusive, will cause the program to stop 
and an appropriate warning message to appear. 

To insure symmetry of solutions for panels which are symmetric with re- 
gard to stringer locations about y = L /2, care should be taken wtih the input 

«y 

data to see that the plate nodes associated with the stringers are symmetric 
about y = Iy/2. Thus, the number of primary structure finite elements in the 
y direction should not be odd if the number of panel stringers is even. 
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C . INPUT INSTRUCTIONS 


A description of the card input for the IBM 370 and CDC 6600 versions of this 
program is presented in this section. 

In addition to the first two input cards which contain literal data, such as 
special program title and date, in columns 1 through 80 , inclusive; there are six 
groups of input cards containing the following information: 

• Group A - Instructions regarding the type of problem being performed, num- 
ber of iteratious desired, and type of output information- 

• Group B - Details of the geometric configuration and finite element mesh of 
the primary structure and tiles. (Card B.4 is omitted if there are no 
tiles) 

• Group C - Defines the primary structure and RSI temperature dependent 
material properties. If there is no TPS, cards C.3 through C'.ll are omitted. 

• Group D - Specifies the primary structure boundary conditions 

• Group E - Describes the mechanical loading upon the primary structure as well 
as its temperature. These cards are omitted when the vibration option is 
used 

m Group F - Defines the RSI temperature distribution. These cards are 
omitted if there is no TPS. 



B— 5 



PROGRAM OPTIONS AND CONTROL - Sheet 1 of 2 
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CARD(S) COL(S) FORMAT SYMBOLS UNITS DESCRIPTION 



If there are no tiles then n and n , together with n and n , are still required since they determine 

x y B*~ -Ucl 

the pri ma ry structure finite element grid. In analyzing panels without tiles, leave out cards B. 4, C.3 
through C.10 and all "F" cards. 
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Figure B-3 RSI Tile Parameters 



MATERIAL PROPERTIES - Sheet 1 of 2 (See Figure B-4) 
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Figure B-4 Orthotropic Stress-Strain Law for 3-Dimensional Elements 
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TEMPERATURE DEPENDENT MATERIAL PROPERTIED , Sheet 1 of 2 
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TEMPERATURE DEPENDENT MATERIAL PROPERTIES - Sheet 2 of 2 
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BOUNDARY CONDITIONS - Sheet 2 of 3 (See Figure B-5) 
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(IN-PLANE Y DEFLECTION) 
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Figure B-5 Primary Stucture Boundary Condition Notation 



BOUNDARY CONDITIONS - Sheet 1 of 3 (See Figure B-5) 



of plane unit rotation 



BOUNDARY CONDITIONS - Sheet 1 of 3 (See Figure B-5) 
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PRIMARY STRUCTURE LOADING (See Figure B-6) 
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Note: Leave out cards E.l and E.2 if vibration option is used 
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Figure B-6 Possible Static Mechanical Loadings Upon Panel 
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D. DESCRIPTION OF OUTPUT 

Output from a typical run of the RESIST computer program is 

explained below in outline form. References in parentheses refer to 

pages in this Appendix. 

1. Program title and date indicating latest update of program 
version which was run. 

INPUT INFORMATION 

2. Listing of input cards , the first two of which are the title 
assigned to any given run by the user. 

3. User selected input options are listed. 

4. Plate, stringer and tile geometry and specification of finite 
element grids for primary structure and tiles (pp. B.8 - B.l]}. 

5. Plate, stringer, strain isolator and arrestor material properties 
(p. B.12 ). Note, if there is no strain arrestor, RSI or isolator 
material properties may be used for the arrestor. If this is done, 
the thickness dimension of the usual isolator or RSI should be 
appropriately reduced to compensate for this addition. 

6. Temperature -dependent RSI material property data used for generating 
curves used internally by program to compute RSI average finite 
element properties. 

7. Plate and stringer boundary conditions (pp. B.16-B.21). 

8. Applied primary structure static mechanical and thermal loading 
if not a vibration problem. 

9. RSI temperature distribution input data. Used for property data 
(item 6 above) and thermal loading if a statics problem. 

OUTPUT INFORMATION 

1 

10. Map showing typical tiles three dimensional finite element ordering, 
by layers. Top, or first layer also corresponds to two-dimensional 
tile coating elements as well. 
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11 . 


Map showing ordering of a typical tiles finite element nodes 
by layers. 

12. Position and temperatures for a typical tile in a local coordinate 
system (reference Figures 3. a &.b). 

13. Global geometry of primary structure nodes and plate nodal degree - 

of -freedom numbering . I> x , and refer to nodal deflections, 

and R , R and R are the nodal rotations. Nodes with no degrees- 
x’ y z & 

of -freedom are used to define the stringer centroids and axes. 

14. a. Statics Option: Primary structure nodal deflections by iteration 

number. Nodes with the same x coordinate are grouped together. 

These groups are separated with dashed lines. 

14. b. Vibration Option: Mode numbers, approximate frequencies and 

corresponding modal error bound (which should be less than 2$> 
to be a reliable approximate mode). This is followed by the 
primaiy structure mode shapes with a similar nodal deflection for- 
mat as for the Statics Option. 

15 . If requested by the user, the computed convergence parameter 

is printed out along with the input quantity it was tested against. 

This is done for each iteration after the first for the Static 
Option. The primary structure degree-of-freedom with the largest 
change from the previous iteration is also identified. 

16. Tile nodal displacements by tile and iteration number. For a 
vibration option, this calculation and the subsequent ones are 
performed only for the user-specified vibration mode. 

17. Three dimensional tile stresses and strains for the bottom two 
layers of elements by element number. These quantities are 
computed at each element's 8 Gauss integration points. Gauss point 
stresses are believed to be more accurate than nodal values and provide 
more detail than simply the elements average stresses. 
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18. Three dimensional element average stresses and strains (by tile 
and iteration number) . 

19. Two-dimensional element average coating stresses. Coating 
element numbers correspond to three dimensional element numbering 
directly below them. 

20. a. Statics Option: Repeat of items 16-19 for each tile. Repeat 

of item l4.a and 15 for each iteration. 

20. b. Vibrations Option: Computation of Rayleigh Quotient (OMEGA 

SQUARED) if all tiles have been treated. Repeat of items 16-19 
for each iteration. Repeat of items l4.b, 15 and Rayleigh 
Quotient until convergence or last iteration is performed. 

21. Plate element stresses and strains for mid and/or top and/or 
bottom surfaces. This computation is done after each iteration 
if requested by the user. Otherwise, it is computed only after 
convergence or the last iteration is performed. 
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E. SAMPLE PROBLEMS 


Output for three sample problems, one vibration case and 
two statics cases, are presented in the remaining pages of this 
report. Only portions of the output for each problem are 
shown. However, the pages presented are representative of the 
types of information, and their respective formats, which the 
RESIST Program can deliver. 
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OF GRUMMAN AEROSPACE CORPORATION 



PROGRAM LISTING OF INPUT DATA CARDS 
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TEMPER ATUR E DEPENDENT MATERIAL PROPERT I E? 
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EDGE PLATE OUT OF PLANE PLATE IN PLANE STRINGERS 




MO ST AT T C THFPMAL LOADING 
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NODE LOCAL TILE COORDINATES TEMPERATURE 
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APPENDIX C 


FINITE ELEMENT TYPES USED IN RE*S*I*ST 

The active degrees-of-freedom (dof) for an element are defined as the dof for 
which there is at least one nonzero term in the corresponding row of the element 
stiffness matrix. The nodes and active dof per node, for each type element 
included in the RESIST computer program, are defined in (Figures C-l through C-4). 




z 



• BENDING AND MEMBRANE STIFFNESS 

• FOUR NODES AT PLATE MIDDLE SURFACE 

• ACTIVE DOF PER NODE ARE x, y AND z DEFLECTIONS 
AND x AND y ROTATIONS 

• LUMPED MASSF C XT NODES 


Figure C-1 Primary Structure Isotropic Plate Element 
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• ARBITRARY CROSS SECTION BENDING, TWISTING AND 
AXIAL STIFFNESS 

• TWO CENTROIDAL NODES AND TWO ATTACHMENT NODES 

• ACTIVE DOF ARE x, y AND z DEFLECTIONS AND ROTATIONS 
AT ATTACHMENT NODES ONLY. NO ACTIVE DOF AT 
CENTROIDAL NODES 

• LUMPED MASSES AT CENTROIDAL NODES 


Figure C-2 Primary Structure Uniform Stiffener Element 


C— 2 



NODE 



• THREE DIMENSIONAL STIFFNESS ELEMENT 

• EIGHT CORNER NODES 

• ACTIVE DOF PER NODE ARE x, y AND z DEFLECTIONS 

• LUMPED MASSES AT CORNER NODES 


Figure C-3 Three Dimensional Orthotropic TPS Element 



• TWO DIMENSIONAL MEMBRANE ELEMENT 

• FOUR CORNER NODES 

• ACTIVE DOF PER NODE ARE x AND y DEFLECTIONS ONLY 

» INERTIA EFFECTS NEGLECTED_F_OR-THESE-ELEMENTS 


Figure C-4 TPS Thin Membrane Coating Element 
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